#include"potential.h"
#include<gsl/gsl_sf.h>
#include<cmath>
double Malflie_Tjon(double r)
{
  return Vr*exp(-mur*r)/r+Va*exp(-mua*r)/r;
}
double Malflie_Tjon(int l,double k1,double k2)
{
  double xr=(k1*k1+k2*k2+mur*mur)/(2.*k1*k2);
  double xa=(k1*k1+k2*k2+mua*mua)/(2.*k1*k2);
  return (Vr*gsl_sf_legendre_Ql(l,xr) +
	  Va*gsl_sf_legendre_Ql(l,xa) )/(2*k1*k2);
}

complexd Malflie_Tjon(int l,complexd k1,complexd k2)
{
  complexd xr=(k1*k1+k2*k2+mur*mur)/(2.*k1*k2);
  complexd xa=(k1*k1+k2*k2+mua*mua)/(2.*k1*k2);
  complexd x1 = k1 * k1 + k2 * k2 + 2. * k1 * k2 + mur*mur;
  complexd x2 = k1 * k1 + k2 * k2 - 2. * k1 * k2 + mur*mur;
  complexd y1 = k1 * k1 + k2 * k2 + 2. * k1 * k2 + mua*mua;
  complexd y2 = k1 * k1 + k2 * k2 - 2. * k1 * k2 + mua*mua;
  complexd legendre1;
  complexd legendre2;
  if ( l == 0 )
    {
      legendre1 = 0.5 * log(x1/x2);
      legendre2 = 0.5 * log(y1/y2);
    }
  else if( l == 1 )
    {
      legendre1 = ( 0.5 * xr * log(x1/x2) - 1. );
      legendre2 = ( 0.5 * xa * log(y1/y2) - 1. );
    }
  return Vr/( 2. * k1 * k2 )*legendre1 + Va/( 2. * k1 * k2 )*legendre2;
}
